function [V,R,nitr] = Iram(A,which,k,m,f,X_ref,L_ref)
%
%
%   Input:  A     -- an n by n matrix
%
%
%           which -- a string indicating which eigenvalues are
%                      wanted.  Options are:
%
%                         'LR', 'SR', 'LM','LA','SA'
%
%           k    -- a positive integer
%
%           m    -- a positive integer  k < m << n
%
%                     Recommended: m > 2k
%
%           f    -- an nonzero n vector
%
%                     If f does not appear, f <- randn(n,1) will occur
%
%
%   Output: V    -- an n by k orthogonal matrix
%
%           R    -- a k by k quasi upper triangular matrix
%
%           ritz -- a vector containing Ritz estimates for the
%                   k Ritz values (eigenvalues of R);
%
%
%           with   norm(AV - VR)  ~  norm(ritz)
%
%           Hence, a partial real Schur form of order k
%
%
%   D.C. Sorensen
%   2 March 2000
%
[n,n] = size(A);
if (nargin < 5), f = randn(n,1); end

V = zeros(n,m); H = zeros(m,m);

tol = 1e-5; ritz = 1; nitr = 0;

ko = 1;
X_err = 1;
L_err = 1;
msave = m; ksave = k;  
%while (norm(ritz) > 1e-5 & nitr < 400)
while( (X_err > tol || L_err > tol) && nitr < 200)
    m = msave + 4; k = ksave + 0;
    nitr = nitr+1;
    [V,H,f] = Arnold(A,V,H,f,ko,m); 
    [mu, ritz]= select_shifts(H,which);
    Q = eye(m);
    j = m;
    while(j > k)
        %
        %            Apply unwanted eigenvalues as shifts with QR steps
        %
        [Q,H] = qrstep(Q,H,mu(j),1,m);
        %
        %            Decrease j by 2 if m_j imaginary and by 1 if real
        %
        if (abs(imag(mu(j))) > 0),
            j = j-2;
        else
            j = j-1;
        end
    end
    ko = j;
    
    ritz = norm(f)*ritz(1:ksave);
    f = V*Q(:,ko+1)*H(ko+1,ko) + f*Q(m,ko);
    V(:,1:ko) = V*Q(:,1:ko);
    
    % Convergence 
    [QQ,RR] = schur(H(1:ko,1:ko));
    VV      = V(:,1:ko)*QQ;
    et      = diag(RR);             % put eigval's in a vector
    [e,idx] = sort(et,1,'descend'); % sort them, high to low, and save index
    X       = VV(:,idx(1));
    L       = e(1);
    X_err   = norm(X - sign(X(1)/X_ref(1))*X_ref, 2);
    L_err   = abs(L - L_ref);
    if ( mod(nitr, 2) == 0 )
       disp(sprintf(...
        ' Iter = %5i, X_err = %6.5e, L = %6.4e, L_err = %6.3e', ...
        nitr, X_err,L,L_err))
    end
end
%
%  Change basis to partial Schur form
%

[Q,R] = schur(H(1:ko,1:ko));
V = V(:,1:ko)*Q;
